c Resampling to generate std errors for estimate of type proportions.
c
      program rsbay51
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,nrep=199)
c
      character*24 fdate
      dimension px(ndec,2),pxm(ngrid,ngrid,ndec,2)
      dimension w0(ngrid),w(ngrid,ngrid),cum(nplay)
      dimension gi(nplay,ngrid,ngrid),gt(ngrid,ngrid)
      dimension f0(nplay,ngrid,ngrid),g(ngrid,ngrid)
      dimension cat10(ngrid,ngrid),q10(4)
      dimension b10(4),b20(4),z10(4),z20(4)
      dimension is0(nplay,ndec),is(nplay,ndec)
      common /mats/is0,is
      common /pmatrx/pxm
c
      open(unit=11,file='ra2011.data')   !with card sleeves on boards
      open(unit=15,file='cat10.data')
c
      read(11,*)((is0(i,n),n=1,ndec),i=1,nplay)
      read(15,*)((cat10(j,k),k=1,ngrid),j=1,ngrid)
c
c-----------------------------------------------------------------
      write(*,'('' ***** rsbay51 ******'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/" All 2011 Sessions"//)')
c
      do 5 i=1,ngrid
         if ((i.eq.1).or.(i.eq.ngrid)) then
            w0(i) = 0.5d0
         else
            w0(i) = 1.d0
         endif
    5 continue
c
      do 8 id=1,nplay
         do 7 k=1,ndec
            is(id,k) = is0(id,k)
    7    continue
    8 continue
c
      dz=1.d0/float(ngrid-1)
      do 20 k=1,ngrid
         beta = dz*(k-1)
         do 15 j=1,ngrid
            if (j.eq.1) then
               pz = 0.999d0
            else
               pz = 1.d0 - 0.5d0*dz*(j-1)
            endif
            gam = 2.d0*dlog(pz/(1.d0-pz))
            call fpx(gam,beta,px)
            do 10 m=1,ndec
               pxm(j,k,m,1) = px(m,1)
               pxm(j,k,m,2) = px(m,2)
   10       continue
            w(j,k) = w0(j)*w0(k)
   15    continue
   20 continue
c
      dseed = 90001.d0
      irep=0
   25 continue
      do 28 j=1,ngrid
         do 27 k=1,ngrid
            g(j,k) = 0.d0
            gt(j,k) = 0.d0
   27    continue
   28 continue
      do 50 id=1,nplay
         cum(id) = 0.d0
         do 35 j=1,ngrid
            do 30 k=1,ngrid
               call fam(id,j,k,t)
               f0(id,j,k) = t
               cum(id) = cum(id) + w(j,k)*t
   30       continue
   35    continue
         do 45 j=1,ngrid
            do 40 k=1,ngrid
               gt(j,k) = gt(j,k) + f0(id,j,k)/cum(id)
   40       continue
   45    continue
   50 continue
c
      xn = 1.d0/float(nplay)
      tt=0.d0
      qt10=0.d0
      do 55 i=1,4
         q10(i)=0.d0
         b10(i)=0.d0
         z10(i)=0.d0
   55 continue
      do 80 id=1,nplay
         cumi = 0.d0
         do 65 j=1,ngrid
            do 60 k=1,ngrid
               pjk = gt(j,k) - f0(id,j,k)/cum(id)
               gi(id,j,k) = f0(id,j,k)*pjk
               cumi = cumi + w(j,k)*gi(id,j,k)
   60       continue
   65    continue
         do 75 j=1,ngrid
            do 70 k=1,ngrid
               g0 = xn*w(j,k)*gi(id,j,k)/cumi
               g(j,k) = g(j,k) + g0
               tt = tt + g0
   70       continue
   75    continue
   80 continue
c
      do 90 j=1,ngrid
         pz = 1.d0 - 0.5d0*dz*(j-1)
         do 85 k=1,ngrid
            beta = dz*(k-1)
            k10 = cat10(j,k)+1
            q10(k10) = q10(k10) + g(j,k)
            b10(k10) = b10(k10) + g(j,k)*beta
            z10(k10) = z10(k10) + g(j,k)*pz
   85    continue
   90 continue
      qt10 = q10(1)+q10(2)+q10(3)+q10(4)
      do 95 i=1,4
         b10(i) = b10(i)/q10(i)
         z10(i) = z10(i)/q10(i)
   95 continue
      write(*,'(i4,5f8.4)')irep,(q10(j),j=1,4),qt10
c      write(*,'(i4,8f8.4)')irep,(b10(j),j=1,4),(z10(j),j=1,4)
      if (irep .eq. nrep) goto 100
c
c begin iterations
      irep=irep+1
      call data(irep,dseed)
      goto 25
  100 continue
c
      end
c***********************************************************************
      subroutine fam(id,j,k,t)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2)
      dimension is0(nplay,ndec),is(nplay,ndec)
      common /mats/is0,is
      common /pmatrx/pxm
c
      t=1.d0
      do 10 n=1,ndec
         m = is(id,n)
         t = t*pxm(j,k,n,m)
   10 continue
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)
c
c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end                                         
c***************************************************************
      subroutine data(it,dseed)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=11,mp=ngrid*ngrid)
      dimension is0(nplay,ndec),is(nplay,ndec)
      dimension pxm(mp,ndec,2)
      common /mats/is0,is
      common /pmatrx/pxm
c
      do 50 id=1,nplay
         r = ggubfs(dseed)
         jd = nplay*r + 1
         do 10 k=1,ndec
            is(id,k) = is0(jd,k)
   10    continue
c
c      write(20,'(i3,x,i2,x,i2,x12i3)')it,id,jd,(is(id,k),k=1,ndec)
   50 continue
      return
      end
c
c========================================================================
      double precision function ggubfs (dseed)
      double precision   dseed
      double precision   d2p31m,d2p31
      data              d2p31m/2147483647.d0/
      data              d2p31 /2147483648.d0/
      dseed = dmod(16807.d0*dseed,d2p31m)
      ggubfs = dseed / d2p31
      return
      end



